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I.    INTRODUCTION 

A.    BACKGROUND 

A  primary  goal  of  a  computer  graphics  system  is  to  provide  the  user  with 
different  views  of  objects.  Sometimes  the  objects  that  the  user  manipulates  are 
simple  in  nature  and  can  be  constructed  easily  with  the  primitives  provided  by 
the  graphics  support  package.  However,  in  most  applications  areas  the  user  is 
concerned  with  more  complex  objects.  The  display  of  three-dimensional  surfaces  is 
one  such  application  area  in  computer  graphics  and  is  the  area  that  this  study 
explores. 

It  is  the  primary  intent  of  this  study  to  stimulate  the  reader's  interest  in  the 
area  of  three-dimensional  surface  generation  and  display.  To  provide  this 
stimulation,  we  combine  the  power  of  certain  mathematical  techniques  and  a  high 
performance  graphics  environment  to  design  and  implement  a  set  of  functions 
that  can  be  used  to  create,  manipulate,  and  display  three-dimensional  solid-filled 
surfaces.  Once  developed,  the  reader  will  not  only  be  able  to  use  these  functions 
to  explore  the  design,  representation,  and  rendering  of  such  surfaces  but  also  will 
be  able  to  use  these  functions  in  other  fields  that  can  derive  benefit  from  their  use 
such  as  cartography,  robotics,  computer  vision,  and  artificial  intelligence. 


B.    PROJECT  DESCRIPTION 

1.  Motivation 

For  developing  graphics  applications,  the  typical  graphics  environment 
provides  the  user  with  a  variety  of  sophisticated,  powerful  tools  and  a  high-level 
language.  Almost  inevitably,  an  application  calls  for  something  that  is  not 
provided  for  or  if  provided  is  not  acceptable.  The  lack  of  solid-filled  curved 
surface  support  in  one  particular  high-performance  graphics  workstation 
prompted  this  study.  By  providing  this  capability  we  can  create  applications  that 
enhance  and  expand  what  we  are  able  to  do  and  conceive  of  doing  on  these 
workstations. 

2.  Proposed  Capabilities 

To  provide  this  capability  for  generating  solid-filled  parametric  bicubic 

surface     patches,    we     must     look     for     a    method     that     is    simple,    powerful, 

understandable,  and  implementable  in  software.    One  approach  would  be  to  start 

from  scratch  and  develop  a  solution.    This,  however,  is  usually  not  a  wise  way  to 

proceed.    A  solution  derived  in  this  manner  is  likely  to  have  severe  limitations 

such  as  being  computationally  inefficient,  lacking  robustness,  and  difficult  to  use. 

Another   more    prudent    way    to    proceed    is    to    find    some    basic    mathematical 

techniques  that  can  be  applied  to  the  problem  in  a  form  consistent  with  the 

existing   graphics   environment.     Solutions   developed    in   this   manner  are  more 

easily   accepted   since   their   basis   lies   in   proven   mathematical    techniques   and 

generally     meets     the     requirements     we     seek     of     being     simple,     powerful, 
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understandable   and    implementable.     We   chose  to   develop  our  new   capability 
following  the  latter  method  -  using  pre-existing  mathematics. 

a.  Parametric  Bicubic  Surface  Construction 

One  of  the  primary  capabilities  needed  when  working  with  surfaces  is 
to  have  an  effective  method  for  describing,  manipulating  and  displaying  them. 
For  this,  we  base  our  work  with  surfaces  on  a  mathematical  model  -  the 
parametric  bicubic  surface  patch.  Choice  of  this  method  allows  the  user  to  rapidly 
develop  a  smooth  surface  and  display  it  as  a  wireframe  image  through  the 
specification  of  only  a  few  points  called  control  points. 

b.  Polygonal  Parametric  Bicubic  Surface  Decomposition 

Another  capability  required  is  to  extract,  from  the  mathematical 
model,  the  information  needed  to  construct  a  solid  surface.  For  this,  we  need  to  be 
able  to  decompose  any  surface  represented  in  parametric  bicubic  form  into 
arbitrarily  small  polygons.  These  polygons  can  then  be  used  in  the  construction 
of  a  polygon  mesh.  This  capability  permits  the  user  to  use  standard  or 
customized  algorithms  to  manipulate  the  surface. 

C.  PROGRAMMING  ENVIRONMENT 

The  IRIS  Turbo  2400  Graphics  Workstation,  manufactured  by  Silicon 
Graphics,  Inc.,  is  the  target  programming  environment  for  the  design, 
development,  and  implementation  of  the  parametric  bicubic  surface  constructor 
and   polygonal   bicubic   surface   decomposition  functions.   The   IRIS   has  special- 
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purpose  hardware  that  is  designed  to  replace  less  efficient  software.  The  system 
supports  real-time  color  graphics,  the  Unix  operating  system  software,  and  the  C 
programming  language.  A  high  resolution  color  monitor  provides  the  output 
device  for  displaying  the  graphical  output  of  the  modeling  and  user  defined 
functions. 

In  addition  to  the  standard  programming  support  provided  by  the  UNIX 
operating  system,  the  IRIS  has  an  extensive  collection  of  utility  and  graphics 
functions  contained  in  a  Graphics  Library.  This  library  provides  the  user  high- 
level  access  to  the  hardware,  enabling  graphical  objects  to  be  easily  manipulated 
as  geometrical  objects  (points,  lines,  polygons,  etc.)  rather  than  pixels.  A  series  of 
coordinate  systems  and  mapping  instructions  also  provides  the  user  with  the 
capability  to  define  such  objects  in  world  coordinate  space. 
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II.   THREE  DIMENSIONAL  REPRESENTATION  OF  SURFACES 

We  stated  in  the  previous  chapter  that  an  important  application  area  in 
computer  graphics  is  concerned  with  the  three-dimensional  representation  of 
surfaces  but  we  did  not  describe  how  one  might  do  this.  There  are  many  ways  to 
represent  surfaces.  The  two  most  commonly  used  representations  are:  polygon 
meshes  and  parametric  bicubic  patches. 

A.    POLYGON  MESHES 

A  polygon  mesh  is  nothing  more  than  a  set  of  connected  polygonally  planar 
surfaces.  There  are  several  ways  in  which  these  planar  polygons  can  be 
represented  -  explicit  polygons,  vertex  list  pointers,  explicit  edges,  etc..  Each  of 
these  representations  has  its  own  advantages  and  disadvantages  and  various 
criteria  can  be  applied  to  evaluate  the  representation.  Criteria  such  as  how  much 
primary  and  secondary  storage  is  available,  how  easy  is  it  to  identify  the  polygons 
sharing  an  edge,  how  difficult  is  it  to  display  the  mesh,  name  only  a  few 
commonly  used  for  evaluation.  Regardless  of  how  the  mesh  is  represented,  there 
are  many  algorithms  available  for  processing  it.  For  example,  algorithms  have 
been  developed  to  remove  hidden  surfaces  from  an  object,  while  others  can 
produce  lighting  and  shading  effects  on  the  surface  of  the  object. 
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Many  objects  that  have  planar  features  such  as  buildings,  tables,  and  cabinets 
can  easily  be  represented  by  a  polygon  mesh.  However,  polygon  meshes  are  not 
limited  to  representing  only  planar  objects.  They  can  also  be  used  to  represent 
curved  surfaces.  To  represent  a  curved  surface  with  a  polygon  mesh,  one  creates 
an  approximation  to  the  surface  by  using  arbitrarily  small  polygons.  The  smaller 
the  polygons  the  better  the  approximation.  However,  certain  difficulties  arise 
when  representing  surfaces  as  polygon  meshes.  As  one  approximates  the  surface 
with  smaller  and  smaller  polygons,  both  space  and  execution  time  of  the 
algorithms  that  process  the  mesh  increase  linearly. 

B.    PARAMETRIC  BICUBIC  PATCHES 

Parametric  bicubic  patches  are  mathematical  models  of  a  surface  and 
represent  one  of  the  simplest  mathematical  elements  we  can  use  to  model  an 
arbitrary  surface.  A  patch  can  be  defined  as  a  curve  bounded  collection  of  points 
whose  coordinates  are  given  by  a  continuous,  two-parameter,  single  valued 
function  in  the  form 

x  =  x(s,t) 
y  =  y(s,t) 

z  =  z(s,t) 
where  the  parameters  are  restricted  to 

s,t  e  [0,11. 
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Patches  are  used  to  describe  the  surface  of  an  object  in  a  piecewise  manner 
much  like  the  polygon  mesh.  That  is,  many  patches  are  used  to  describe  the 
surface  of  an  object  and  the  total  surface  is  generated  by  displaying  all  of  its 
individual  patches. 

Bicubic  patches  are  most  often  used  to  generate  line  drawings  of  three- 
dimensional  objects  -  often  called  wireframe  representations.  One  benefit  of  using 
bicubic  patches  is  that  fewer  bicubic  patches  than  polygonal  patches  are  needed 
to  represent  a  curved  surface  to  a  given  accuracy.  On  the  other  hand,  the 
algorithms  for  working  with  bicubic  surfaces  are  more  complex  than  those  for 
polygon  meshes.  In  addition,  wireframe  representations  can  exhibit  certain 
deficiencies.  For  example,  objects  that  are  represented  by  wireframes  are  often 
ambiguous.  That  is,  you  can  look  at  the  same  object  and  get  different  visual 
interpretations.  One  reason  for  the  different  interpretations  is  the  ability  to  see 
through  a  wireframe  representation  of  an  object.  While  not  a  problem  for  some 
applications,  it  can  be  a  serious  problem  for  others. 
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III.    PARAMETRIC  CUBIC  CURVES  AND  SURFACES 

A.    GENERAL 

The  method  that  we  use  to  generate  surfaces  is  based  on  using  parametric 
cubic  curves  so  it  is  helpful  to  review  the  mathematical  basis  for  these  curves.  A 
parametric  cubic  curve  has  the  property  that  x,  y,  and  z  can  be  defined  as  third- 
order  polynomials  for  some  variable  t  [Ref.  l:pp.  34]  : 

x(t)  =  axts  +  bxt2  +  cxt  +  dx, 
y(t)  =  ayts  +  byt2  +  cyt  +  dy, 
z(t)  =  arts  +  brt2  +  crt  +  dt. 

These  equations  are  known  as  the  algebraic  form  of  a  parametric  cubic  curve. 
In  this  form,  we  can  identify  a  unique  set  of  12  constant  algebraic  coefficients. 
These  algebraic  coefficients  determine  a  unique  parametric  cubic  curve;  they 
determine  the  size  and  shape  of  the  curve  and  its  position  in  three-dimensional 
space.  Two  curves  of  the  same  shape  have  different  algebraic  coefficients  if  they 
occupy  different  positions  in  three-dimensional  space.  Because  we  want  to  deal 
with  a  finite  segments  of  the  curve,  we  limit  the  range  of  the  parameter,  without 
loss  of  generality,  to  O^t^l.  We  call  these  finite  pieces  curve  segments.  A  curve 
segment  is  nothing  more  than  a  point-bounded  collection  of  points.  In  our  case 
the  points  are  bounded  at  t=0  and  t  =  l.  Each  equation  in  the  algebraic  form  can 
be  expressed  as  a  vector  product  as  follows: 
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x(t)  =  [t8  t2  t  l]    ^x      =  [t8  t2  t  l]  Cx  =  TCX. 
[dx 

Using  the  vector  product  form  is  usually  more  convenient  as  it  separates  the 
distinct  parameters  of  the  parametric  equation  into  unknown  coefficients  of  x(t) 
and  the  parameter  t  that  we  wish  to  manipulate.  Here  T  is  the  row  vector  of 
powers  of  t,  while  Cx  is  the  column  vector  of  coefficients  of  x(t).  Similarly  the 
parametric  equations  for  y(t)  and  z(t)  can  be  written  as  y(t)  =  TCy  and  z(t)  =  TCt. 
By  varying  the  parameter  t  from  0  to  1  in  each  equation  we  define  the  curve 
segment. 

Arbitrarily  assigning  values  to  these  unknown  coefficients  results  in  defining  a 
curve  in  three-dimensional  space.  However,  it  is  not  easy  to  determine  the 
properties  of  this  curve.  What  we  wish  to  do  is  establish  some  constraints  on 
these  coefficients.  We  want  the  curves  we  generate  to  have  some  predictable 
properties.  To  solve  the  equations  for  these  unknown  algebraic  coefficients,  we 
establish  a  set  of  constraints,  thereby  defining  a  unique  cubic  curve  with 
predictable  properties.  To  illustrate  this  process  we  look  at  some  example  cubic 
curves  for  which  the  constraints  are  well-known. 

B.    CUBIC  CURVE  EXAMPLES 

1.     Hermite  Curve 

The  Hermite  cubic  curve  is  determined  from  its  endpoints  (Pt,  P2)  and 

endpoint  tangents  (R1?  R2).    In  the  literature  [Ref.  2:pp.  516-519]  [Ref.  3:pp.  123- 
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129],  we  find  the  geometric  form  of  the  Hermite  curve  to  be 

Qh(t)  =  TMhGh. 

In  geometric  form.  T  is  the  row  vector  of  the  powers  of  the  parametric 
variable  t,  Mh  is  the  basis  matrix,  and  Gh  is  the  geometry  vector.  A  basis  matrix 
refers  to  a  constraint  procedure  that  is  embodied  in  matrix  form  and  a  geometry 
vector  contains  the  control  points  used  to  guide  the  curve. 

In  this  particular  case,  the  Hermite  basis  matrix  (Mh)  is 


2-211 

-3      3     -2    -1 

0       0       10 

10       0       0 


and  the  Hermite  geometry  vector  is 


Pi 
P2 

Ri 

if  •> 


Now  using  the  above  formulation,  given  two  points  and  their  tangents,  we 
can  evaluate  x(t),  y(t),  and  z(t)  for  O^t^l  and  find  all  points  on  the  Hermite  form 
of  the  cubic  curve  from  Px  to  P2  with  starting  tangent  vector  Rx  and  ending 
tangent  vector  R2.  It  is  through  these  constraints  (Mh)  that  the  control  points 
(Gh)  control  the  parametric  equations  and  produce  an  equation  that  can  generate 
a  discretely  sampled  curve  segment  in  three-dimensional  space. 
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If  we  take  the  product  TMh,  we  have 

TMh  =   [(2t8-3t2  +  l)  (-2ts  +  3t2)  (t8-2t2  +  t)  (t5-t2)  J. 

These  four  functions  of  t  in  the  product  TMh  are  often  called  blending  functions 
[Ref  l:pp.  48-52].    As  the  name  implies,  they  blend  the  effects  or  contributions  of 
the  endpoints  and  tangent  vectors  to  produce  the  intermediate  point  coordinate 
values  over  the  domain  of  t. 
2.     Bezier  Curve 

The  defining  form  for  a  Bezier  cubic  curve  is  similar  to  the  Hermite  form. 
The  difference  is  in  the  definition  of  the  endpoint  tangent  vectors.  The  Bezier 
form  uses  four  points  (  Pl,  P2.  P8,  P4  )  instead  of  2  points  and  2  tangent  vectors. 
The  tangent  vectors  at  the  endpoints  in  Bezier  form  are  determined  by  the  line 
segments  PjP2  and  PSP4  •  The  Bezier  cubic  curve  passes  through  the  first  and 
fourth  control  points  (Pj  and  P4)  and  uses  the  second  and  third  points  (P2  and 
P,)  to  determine  the  shape  of  the  curve.  [Ref.  l:pp.  113-125]  [Ref.  2:pp.  519-521]. 

The  geometric  form  of  the  Bezier  curve  is 

Qb(t)  =  TMbGb 
where  the  Bezier  basis  matrix  (Mb)  is 


My 


-1  3-3  1 
3-630 

-3300 
10       0      0 
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and  the  Bezier  geometry  vector  is 


P, 
P2 
P8 
P< 


The  Bezier  form  of  the  cubic  curve  is  more  widely  used  in  computer 
graphics  than  the  Hermite  form.  A  primary  reason  for  its  popularity  is  that  the 
geometry  matrix  of  four  points  (Gb)  is  more  intuitive  for  an  interactive  user.  The 
user  has  only  to  manipulate  the  four  points  and  does  not  have  to  specify  the 
tangent  vectors.  It  is  usually  easier  for  a  person  to  think  about  manipulating 
points  rather  than  trying  to  manipulate  points  and  tangent  vectors. 
3.     Other  Useful  Cubic  Curves 

The  Hermite  and   Bezier  cubic   curves  are  not  the  only  forms  of  cubic 
curve  that  are  available.    Two  others  are  the  Cardinal  Spline  and  the  B-Spline. 

a.     Cardinal  Spline 

The  Cardinal  Spline  curve  passes  through  the  two  interior  control 
points  (P2  and  P8)  and  uses  the  points  Pj  and  P4  to  define  the  shape  of  the  curve 
[Ref.  4:p.  11-4].    The  geometric  form  of  the  Cardinal  curve  is 

Qe(t)  =  TMcGe 

where  the  Cardinal  basis  matrix  (Mc)  is 


18 


—a 

2-a      a-2 

a 

2a 

a-3    3-2a 

—a 

-a 

0           a 

0 

0 

1            0 

0 

and  the  Cardinal  geometry  vector  is 


P, 
P2 
Ps 
P4 


The    scalar   coefficient    a    in    the    Cardinal    basis    matrix   must    be   positive   and 
determines  the  length  of  the  tangent  vector  at  point  P2  and  Ps. 
b.     B-Spline 

The  geometric  form  of  the  B-Spline  curve  is 


Qbi(t)  =  TMb,G 


bs 


where  the  B-Spline  basis  matrix  (Mb8)  is 


-1      3-3     1 

3-630 

-3030 

14       10 


and  the  B-Spline  geometry  vector  is 


Pi 
P2 
P8 
P4 
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In  general,  the  B-Spline  curve  does  not  pass  through  any  control 
points  but  is  continuous  and  also  has  continuity  of  tangent  vectors  and  of 
curvature  (that  is,  first  and  second  derivatives  are  continuous  at  the  endpoints). 
The  Hermite  and  Bezier  forms  have  only  first-derivative  continuity  at  the 
endpoints,  but  do  pass  through  control  points  [Ref.  l:pp.  125-146]  [Ref.  2:pp. 
521-523]. 

C.    DEFINING  SURFACES 

By  adding  a  new  parameter  s  and  additional  algebraic  coefficients  to  the  cubic 
curves  in  the  previous  section,  we  can  define  the  algebraic  form  of  a  bicubic 
surface  patch  [Ref.  l:pp.  156]  as 

j=s j=s 

P(s,t)  =  £  SauS'tJ 
1=0 j=o 

with  the  restriction  on  the  parametric  variables  to 

S,t£    [0,1]- 

By  varying  both  parameters  from  0  to  1  in  each  equation,  we  define  all 
points  on  the  surface  patch.  Assigning  one  parameter  a  constant  value  and 
varying  the  other,  results  in  a  cubic  curve. 

Expanding  the  above  equation  in  terms  of  x(s,t)  and  noting  that  the  terms  for 
y(s,t)  and  z(s,t)  are  similar  we  have 
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x(s.t)  =  aS8sst8  +  a82sst2  +  a81s8t  +  asoss  +  a3Ss2ts  +  a22s2t2  +  a21s2t  +  a20s2 
+  alssts  +  a12st2  +  anst  +  a10s  +  aosts  +  a02t2  +  a01t  +  a00. 


Written  in  vector  product  form 


x(s,t)  =  SC^T1 


where 


S  =  [ss  s2  s  1  ]  , 
T  =  It8  t2  t  1  I, 


Cv  = 


aSS  a82  a81  a80 

a2S  a22  a21  a20 

a18  a12  an  aio 

a08  a02  a01  a00 


and  TT  is  the  transpose  of  the  matrix  T  . 

From  these  equations  we  can  see  that  there  are  48  degrees  of  freedom  or 
algebraic  coefficients  that  we  must  specify.  Like  the  cubic  curve,  a  change  in  any 
one  of  these  coefficients  defines  a  different  surface. 

The  complete  algebraic  manipulation  of  the  equations  to  arrive  at  the 
following  equation  is  similar  to  that  of  the  curve  process  described  in  the  previous 
section.  For  the  Bezier  surface  patch,  the  geometric  form  of  the  equation  is: 


x(s,t)  =  SMbQxMbTT 


TtT 


where  Mb  is  the  same  Mb  as  in  the  Bezier  curve  equation,  Mb  is  its  transpose,  and 
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Qx   is  the  x  component  of  the  sixteen  control  points  of  a  surface  patch.     The 
matrix  Qx  is 


Qx  = 


Pix  Psx  P9x  P18x 

P2,  P6X  P10x  P14x 

P«x  Prx  Pllx  P15, 

P<.  P8,  Pl2,  Pl6, 


and  similarly  for  Qy  and  Qt. 

Since  we  must  provide  three  4x4  matrices,  one  for  each  of  component  x,  y, 
and  z,  it  can  be  seen  that  we  have  specified  the  48  degrees  of  freedom  as  in  the 
algebraic  form. 

As  can  be  seen  by  the  above  equations,  a  bicubic  surface  patch  can  be  defined 
by  a  set  of  16  control  points  and  a  basis  matrix.  By  manipulating  the  control 
points,  we  change  the  shape  of  the  surface  as  constrained  by  the  basis  matrix. 
We  take  this  knowledge  with  us  as  we  design  our  functions  in  the  next  chapter. 
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IV.   DESIGN  AND  IMPLEMENTATION  OF  OUR  SURFACE  FUNCTIONS 

A.    OVERVIEW 

Having  the  mathematics  developed  in  the  previous  chapter  does  us  no  good  if 
we  can  not  to  put  it  to  use.  In  our  case,  this  means  being  able  to  draw  a 
parametric  bicubic  surface.  It  is  at  this  point  that  we  begin  to  see  how  the 
mathematics  can  be  combined  with  the  power  of  the  computer  and  the  graphics 
workstation. 

Up  to  this  point,  we  have  dealt  with  two  forms  of  the  parametric  cubic  curve 
and  parametric  bicubic  surface  -  the  algebraic  form  and  the  geometric  form.  The 
question  now  is  which  one  shall  we  work  with? 

Deciding  what  form  to  use  depends  largely  on  the  application.  If  we  are  given 
or  know  the  algebraic  equations  of  the  curve,  then  the  reasonable  choice  is  the 
algebraic  form.  If  we  plan  to  do  surface  fitting  of  data  or  interactive  design  the 
choice  is  the  geometric.  We  choose  to  use  the  geometric  form.  Our  primary 
reason  for  choosing  it  is  that  the  geometric  form  offers  us  a  greater  insight  into 
the  control  and  behavior  of  curves  and  surfaces  than  is  otherwise  available  with 
the  classical  algebraic  formulation.  It  should  be  noted,  however,  that  it  is 
possible,  through  mathematical  manipulation,  to  convert  from  one  form  to  the 
other  [Ref.  l:pp.  164]. 
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B.  METHODOLOGY 

How  one  proceeds  to  generate  a  surface  impacts  usefulness,  flexibility,  and  the 
ability  to  understand.  As  we  stated  in  the  first  chapter,  the  built-in  functions 
provided  by  a  graphics  environment  are  not  always  exactly  what  we  want.  This 
is  the  case  in  the  IRIS  graphics  environment.  Although  the  IRIS  provides  support 
for  bicubic  surface  patches,  it  does  not  support  surface  patch  decomposition. 
What  we  want  is  the  capability  to  do  both.  We  also  want  this  capability  without 
sacrificing  what  a  user  already  knows  about  how  the  IRIS  supports  bicubic  surface 
patches.  To  achieve  this,  we  have  developed  a  set  of  parallel  routines  [Appendix 
A]  that  provide  nearly  all  the  functionality  as  the  standard  IRIS  functions  while 
at  the  same  time  providing  the  user  with  extended  support  via  three  additional 
functions.  These  extended  functions  allow  the  user  to  have  access  to  the 
triangular  polygons  that  our  new  functions  generate  during  the  construction  of 
the  bicubic  surface  patch.  That  is,  the  IRIS  user  is  able  to  use  the  new  functions 
in  the  same  way  as  he  would  use  the  standard  functions  by  substituting  the 
names  of  the  new  routines  in  place  of  the  standard  IRIS  routines.  If,  however,  the 
user  wishes  to  be  able  to  have  access  to  the  individual  triangular  polygons  that 
make  up  the  surface,  he  has  only  three  additional  routines  to  learn. 

C.  PARALLELING  THE  IRIS  SUPPORTED  FUNCTIONS 

The  IRIS  graphics  environment  has  available  five  functions  for  defining  and 
generating  parametric  bicubic  surface  patches.    Those  five  functions  are: 
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-  dc  [basis  defines  a  basis  matrix 

-  patchbasis  sets  the  current  basis  matrices  for  both  the  s  and  the  t  parametric 
direction 

-  patehcurvcs  sets  the  number  of  curves  used  to  represent  a  patch 

-  patchpreeision  sets  the  precision  at  which  the  curves  are  drawn 

-  patch  draws  the  surface  patch. 

A  complete  description  of  the  functions  and  their  arguments  can  be  found  in 
the  IRIS  Users  Manual  [Ref.  4]. 

To  ease  the  pain  of  learning  new  functions,  our  new  functions  are 
syntactically  identical  to  the  standard  IRIS  functions  with  the  exception  that  the 
new  function  names  are  the  standard  function  names  prefixed  by  the  letter  re. 
These  parallel  functions  are: 

-  ndefbasis 

-  npatehbasis 

-  npatcheurvcs 

-  npatehpreeision 

-  npateh. 

The  usage  and  the  arguments  of  these  parallel  functions  remain  the  same  as 
the  standard  IRIS  functions.  The  only  difference  that  the  user  notices  is  that  the 
wireframe  drawn  looks  like  a  triangular  mesh  instead  of  the  typical  wireframe  and 
that  the  function  npatchprecision  has  no  effect  on  the  displayed  image. 

While  these  routines  seem  to  do  what  the  old  routines  do,  they  are  more 
powerful  because  they  provide  special  extensions  to  the  user.  These  extensions 
provide  the  capability  to  manipulate  the  surface  patch  as  individual  polygons. 
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D.    TRIANGULAR  DECOMPOSITION  OF  SURFACE  PATCHES 

The  extensions  mentioned  in  the  previous  section  are  available  to  the  user  by 
using  three  additional  routines: 

-  Set   User  Routine  for  npateh  provides  an  intercept  function  for  handling  the 
triangular  polygons  generated  in  the  surface  decomposition 

-  User  Routine  provides  the  user  a  way  to  turn  the  intercept  function  on  and 
off 

-  Set  Default  Routine  for  npateh   allows    the   user  to   return   to   the  system 
defined  intercept  function. 

There  is  one  argument  to  the  function  Set  User  Routine  for  npateh.  This 
argument  is  the  name  of  a  user-defined  function  that  expects  to  receive  a  3x3 
array.  This  3x3  array  contains  the  three  vertices  of  a  triangle  where  each  vertex  is 
made  up  of  an  x,  y,  and  z  coordinate.  The  function  User  Routine  expects  one 
argument  also.  If  this  argument  is  zero,  then  the  intercept  function  is  turned  off; 
i.e.,  the  user's  program  cannot  intercept  the  triangles  composing  the  surface. 
Otherwise  the  function  is  activated  allowing  the  user's  program  to  intercept  the 
triangles  comprising  the  surface.  The  function  Set  Default  Routine  for  npateh 
does  not  expect  any  arguments. 

Using  these  functions,  the  user's  program  has  access  to  and  can  manipulate 

the    individual   triangular  components  of  the   surface   patch.     For  example,   an 

individual  surface  patch  can  be  decomposed  into  triangular  polygons  and  then  via 

the  user  intercept  function,  each  polygon  can  be  subjected  to  an  illumination 

model  that  produces  a  realistic  looking  surface  in  three-dimensional  space. 
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The  surface  patch  is  decomposed  into  triangular  polygons  one  at  a  time.  For 
the  user's  program  to  intercept  these  polygons,  the  program  must  have  specified 
an  intercept  function  via  the  Set  User  Routine  for  npateh  and  must  have 
activated  it  via  User  Routine.  Then,  as  each  polygon  is  generated,  the  user's 
function  can  process  them  in  any  way  desired.  They  can  be  stored,  manipulated, 
altered,  etc..  It  is  the  user's  program  that  determines  what  to  do  with  them. 
This  feature  provides  a  tremendous  amount  of  flexibility,  creativity,  and 
applicability  above  what  is  currently  available  in  the  standard  IRIS  graphics 
environment  support  of  surfaces. 

E.    GENERAL  GUIDELINES  FOR  USAGE 

To  prevent  any  unnecessary  problems  in  using  our  new  functions,  we  need  to 
establish  a  basic  set  of  guidelines  or  sequences  of  events  that  should  be  followed. 
If  the  user  does  not  want  to  use  the  special  extensions,  i.e. 
Set  User  Routine  for  npateh  and  User  Routine,  then  a  modified  version  of  the 
standard  IRIS  setup  steps  for  using  surface  patches  can  be  followed.  These  steps 
are: 

-  define  the  appropriate  curve  bases  using  the  ndefbasis  function; 

-  select  a  basis  for  the  s  and  t  parametric  directions  using  the   npatehbasis 

function; 

-  select  the  number  of  curve  segments  to  be  drawn  in  each  parametric  direction 
using  the  npateheurves  function; 

-  draw  the  surface  by  using  the  npateh  function. 

The  only  change  to  the  standard  IRIS  setup  is  that  it  is  not  necessary  to  use 
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the   npatehprceision  function.   This  function  does  not  effect  the  displayed  image 
and  its  only  purpose  is  to  maintain  consistency  with  the  standard  IRIS  functions. 

If  the  user  wishes  to  use  the  extensions  via  the  Set   User  Routine  for  npateh 
and  the  User  Routine  functions,  then  the  following  steps  must  added: 

-  An  intercept  function  must  be  declared  and  defined  before  calling  the 
Set  User  Routine  for  npateh  function.  This  intercept  function  must  be 
declared  as  a  function  returning  an  integer  value  (even  though  it  is  not  used) 
and  must  be  defined  as  receiving  a  3x3  matrix  of  floating  point  numbers, 
where  each  row  contains  one  set  of  x,  y  and  z  coordinates  of  an  intercepted 
triangles  vertex.  The  name  of  this  intercept  function  will  be  the  argument 
given  to  the  Set    User  Routine  for  npateh  function; 

-  Activating/deactivating  the  intercept  function  via  the  User  Routine  function 
can  be  performed  any  time  after  the  above  step  has  been  completed. 

By  using  these  guidelines,  a  user  should  not  have  any  difficulty  in  using  these 

functions.    As  we  will  show  in  the  next  chapter,  these  functions  are  easy  to  use, 

efficient,    and    can    provide    some    impressive    results    when    a    carefully    chosen 

intercept  function  is  used. 
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V.   USAGE  AND  PERFORMANCE 

Having  taken  a  brief  tour  of  the  functions  we  designed  in  the  previous 
chapter,  we  need  to  provide  some  concrete  examples  of  their  usage,  performance 
levels,  and  limitations. 

A.    SAMPLE  PROGRAMS 

Appendix  B  gives  the  listing  for  four  sample  programs  using  the  new 
functions.  Each  program  illustrates  how  the  new  functions  can  be  integrated  into 
the  IRIS  graphics  programming  environment.  Program  #1  draws  2  surface 
patches  in  wireframe  representation.  One  surface  is  drawn  using  the  standard 
IRIS  functions  while  the  other  is  drawn  using  the  parallel  functions.  When  this 
program  is  run  the  user  notices  the  different  appearance  of  the  wireframe  surface 
patch  drawn  with  the  new  functions.  It  has  the  triangular  mesh  appearance 
described  in  the  previous  chapter.  Program  #2  shows  how  a  user-defined 
intercept  function  can  be  used  via  the  three  extension  functions  we  have  designed. 
This  program  intercepts  the  triangles  generated  during  the  patches  decomposition 
and  puts  them  into  an  IRIS  graphical  object.  Program  #3  shows  how  the  user- 
defined  intercept  function  can  be  dynamically  changed  during  execution  and 
Program  #4  shows  how  a  well  chosen  intercept  function  can  be  used  to  produce 
realistic  lighting  of  a  surface. 
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B.    PERFORMANCE  COMPARISONS 

Knowing  that  our  functions  work,  we  would  like  to  know  how  efficient  they 
are.  The  way  that  we  approach  this  question  is  to  compare  our  new  functions  to 
the  standard  IRIS  functions.  Because  we  have  designed  our  functions  to  be 
substitutable  as  a  set  for  the  standard  functions,  we  do  not  have  any  problems  in 
testing  the  relative  performance  of  the  sets.  However,  two  words  of  caution  are  in 
order  before  comparing  these  two  sets  of  functions.  First,  the  IRIS  implements 
many  of  its  graphics  primitives  via  special  purpose  hardware.  This  is  the  case 
with  the  function  patch.  Therefore,  its  parallel  function  npateh,  which  is 
implemented  in  software  is  not  as  fast.  Second,  the  new  function  npatehprceision 
does  not  affect  the  computation  whereas  the  IRIS  function  patehprecision  does 
affect  the  computations.  Keeping  these  points  in  mind,  we  have  developed  a 
simple  benchmark  program,  listed  in  Appendix  C,  that  we  use  to  draw  a 
wireframe  representation  of  a  surface  patch  100  times.  By  executing  this  program 
10  times  and  getting  the  average  times,  we  can  get  an  idea  of  the  performance  of 
the  parallel  set  of  surface  patch  functions  as  compared  to  the  standard  IRIS 
surface  patch  functions. 

The  way  that  we  measure  performance  of  a  particular  program  is  to  use  the 
UNIX  time  command.  The  time  command  returns,  on  program  completion,  the 
time  in  seconds  for  system  time,  user  time,  and  elapsed  time. 

The  benchmark  program  was  executed  in  two  different  modes.    In  the  first 

mode,    the  program  was  executed  without  the  assistance  of  the  IRIS's  floating 
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point  accelerator  (FFP)  while  in  the  second  mode,  the  program  was  executed  with 
the  FFP.  The  results  indicate  that  the  standard  IRIS  functions  are  350%  faster 
without  the  FFP  and  260%  faster  with  the  FFP.  As  expected,  the  new  functions 
are  slower,  however  considering  that  they  were  not  designed  to  replace  the  IRIS 
functions,  these  results  are  good.  Normally,  one  can  expect  an  order  of  magnitude 
increase  in  performance  when  special  hardware  is  used. 

C.    LIMITATIONS 

Nothing  that  can  be  developed  is  without  limitations.  The  reader  should 
recall  that  it  was  certain  limitations  of  the  existing  IRIS  system  that  motivated 
this  study.  The  functions  we  have  designed  and  have  implemented  have  allowed 
us  to  overcome  certain  limitations  in  the  IRIS  graphics  environment.  At  the  same 
time,  these  functions  have  their  own  limitations. 

The  primary  limitation  of  our  parallel  functions  is  speed.  While  these 
functions  have  been  carefully  implemented  using  efficient  algorithms  and  data 
structures,  they  are  not  as  fast  as  using  special  purpose  hardware.  Another 
limitation  deals  with  the  use  of  memory.  The  npatch  function  allocates  memory 
to  save  each  point  on  the  surface  patch.  The  number  of  points  that  are  generated 
are  proportional  to  the  product  of  the  desired  number  of  curve  segments  in  the  s 
and  t  parametric  directions.  For  example,  to  draw  a  surface  patch  with  10  curves 
in  the  s  direction  and  10  curves  in  the  t  directions  requires  at  least  enough 
memory  to  store  300  floating  point  numbers  (one  for  each  x,  y,  and  z  component). 
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To  draw  a  100x100  surface  patch  requires  enough  memory  to  store  30,000  floating 
point  numbers.  Assuming  a  floating  point  number  requires  4  bytes,  the  10x10 
patch  requires  1.17  Kilobytes  of  memory,  while  the  100x100  patch  requires  117.1 
Kilobytes  of  memory. 
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VI.   RECOMMENDATIONS  AND  CONCLUSIONS 

A.    DIRECTIONS  FOR  FURTHER  STUDY 

Bicubic  surface  display  and  generation  is  an  area  of  research  in  computer 
graphics  that  is  exciting  and  important.  Since  the  development  of  high- 
performance  graphics  systems,  the  demand  for  realism  and  real-time  has  increased 
significantly.  Consequently,  the  need  is  great  for  continued  creativity  and 
exploration  in  the  area. 

1.     Development  of  Application  Programs 

The  power  of  these  parallel  surface  functions  we  have  created  can  only  be 
derived  through  the  use  of  the  intercept  functions.  Whether  they  will  be  used  to 
experiment  with  lighting  and  shading  models  or  applied  to  fractal  geometry  can 
only  be  answered  by  time.  However,  it  is  through  creative  experimentation  that 
these  questions  can  be  answered.  Some  areas  for  further  work  are: 

-  Surface- fit  ting  sampled  data 

Surface-fitting  is  the  process  of  constructing  a  representation  to  model  the 
surface  of  an  object  based  on  a  fairly  large  number  of  given  data  points.  By 
taking  these  points  and  chosing  an  appropriate  set  of  surface  constraints,  one 
can  accurately  reconstruct  the  surface.  For  example,  during  this  work,  we 
were  given  a  set  of  digitized  x,  y  and  z  coordinates  for  a  human  head.  By 
successively  extracting  control  points  from  the  data,  we  were  able  to 
reproduce  and  display  the  head  quite  accurately. 
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-  Data  Reduction 

In  many  instances  it  is  possible  to  reduce  the  amount  of  data  needed  to 
properly  reconstruct  a  surface.  For  example,  consider  geographical  terrain. 
Terrain  that  is  relatively  flat  can  be  reconstructed  with  fewer  surface  patches 
that  can  mountainous  terrain.  The  problem  is  that  most  terrain  is  sampled 
at  discrete  intervals,  such  as  every  100  meters,  whether  it  is  flat  or  not.  By 
applying  some  form  of  an  Adaptive  Subdivision  Algorithm  [Ref.  5]  one  can 
reduce  the  amount  of  primary  and  secondary  storage  while  at  the  same  time 
provide  increased  performance  for  display. 

-  Lighting  Models 

Because  the  user  can  intercept  individual  polygons  comprising  the  surface,  it 
is  possible  to  subject  each  polygon  to  a  lighting  model.  While  we  have 
provided  a  simple  example  of  this,  more  sophisticated  lighting  models  could 
be  easily  integrated  through  these  parallel  functions. 

-  Realistic  3-D  Objects 

The  surfaces  of  many  vehicles  such  as  automobiles,  aircraft,  and  ships  can  be 
constructed  with  bicubic  surface  patches.  For  example,  constructing  an 
object  with  surface  patches  and  applying  a  lighting  and  shading  model,  one 
could  develop  an  ship  identification  training  system.  Such  a  training  system 
would  be  a  valuable  asset  in  military  training  environments,  allowing  the 
trainee  to  view  a  particular  class  of  ship  from  any  viewing  angle. 

2.     Improvement  of  Performance 

Real-time    computer    graphics    requires    efficient    algorithms    and    data 

structures.    While  these  functions  were  coded  to  be  as  efficient  as  possible,  while 

preserving    understandability,    there    is    always    room    for    improvement.      One 

suggestion  we  have  is  to  contact  the  developers  of  the  IRIS  graphics  package  for 

insights  into  improving  our  packages  performance.     Such  contact  may  provide 

access  to  low-level  graphic  system  routines  and  techniques  that  could  dramatically 

improve  performance. 
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B.    CONCLUSIONS 

This  study  introduced  the  reader  to  the  world  of  parametric  bicubic  surfaces. 
To  do  this,  we  provided  some  necessary  definitions,  terminology  and  mathematics. 
We  also  designed  and  implemented  a  set  of  software  functions  that  take 
advantage  of  the  information  and  given  them  to  the  reader  for  experimentation. 
The  benefit  that  can  be  derived  from  the  use  of  these  functions  can  only  be 
determined  by  the  passage  of  time. 
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APPENDIX  A     -    FUNCTION  SPECIFICATIONS 


NAME 

ndefbasis  -  defines  a  basis  matrix 

SPECIFICATION 

ndefbasis  (id,  mat) 
long  id; 
Matrix  mat; 

DESCRIPTION 

ndefbasis  allows  the  user  to  define  basis  matrices  for  use  in  the 
generation  of  patches,  matrix  is  saved  and  is  associated  with  id.  id 
may  then  be  used  in  subsequent  calls  to  npatchbasis. 


NAME 

npatchbasis  -  sets  current  basis  matrices 

SPECIFICATION 

npatchbasis (sid,  tid) 
long  sid,  tid; 

DESCRIPTION 


npatchbasis  sets  the  current  basis  matrices  (denned  by  ndefbasis  ) 
for  both  the  s  and  t  parametric  directions  of  a  surface  patch.  The 
current  s  and  t  bases  are  used  when  the  npatch  command  is  issued. 
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NAME 

npatchcurves  -  sets  number  of  curves  used  to  represent  a  patch 

SPECIFICATION 

npatehcurves(scurves,  tcurves) 
long  scurves,  tcurves; 

DESCRIPTION 

npatchcurves   sets  the  current   number  of  s  and   t  curves  used  to 
represent  a  patch  as  a  wireframe. 


NAME 

npatchprecision  -  is  a  null  function. 

SPECIFICATION 

npatchprecision(ssegments,  tsegments) 
long  ssegments,  tsegments; 

DESCRIPTION 

npatchprecision  has  no  functionality  at  the  current  time.  It  is  used 
to  maintain  consistency  with  the  standard  IRIS  function 
patchprecision. 
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NAME 

npatch  -  draws  a  surface  patch 

SPECIFICATION 

npatch(geomx,  geomy,  geomz) 
Matrix  geomx,  geomy,  geomz; 

DESCRIPTION 

npatch  draws  a  surface  patch  using  the  current  npatchbasis  and 
npatchcurves.  The  shape  of  the  patch  is  determined  by  the  control 
points  specified  in  geomx,  geomy,  and  geomz. 


NAME 


Set    User   Routine   for   npatch    -    allows    the    user   to    specify    an 
intercept  function 


SPECIFICATION 


Set    User    Routine   for   npatch(fname) 
int  (*fname)(); 


DESCRIPTION 


Set  User  Routine  for  npatch  allows  the  user  to  set  up  a  function 
that  is  capable  of  intercepting  triangular  polygons  generated  during 
the  decomposition  of  a  surface  patch.  The  number  of  polygons 
generated  is  (  (scurves  -  1)  *  (tcurves  -  1)  *  2.    ). 
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NAME 


User   Routine  -  allows  the  user  to  turn  the  intercept  function  on  and 
off 


SPECIFICATION 


User    Routine  (boolean) 
int  boolean; 


DESCRIPTION 


User  Routine  acts  like  a  switch  allowing  a  user-defined  intercept 
function  to  be  turned  on  and  off.  By  assigning  boolean  the  value  0  the 
intercept  function  is  turned  off.  Integer  values  other  than  0  cause  the 
intercept  function  to  be  turned  on. 


NAME 

Set    Default    Routine    for   npatch  -  resets  the  intercept  function 
to  a  system  defined  default 

SPECIFICATION 

Set    Default    Routine   for   npatch  () 

DESCRIPTION 

Set    Default    Routine   for   npatch  enables  the  user  to  choose  the 
system  defined  intercept  function  poly (3,  Triangle). 
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APPENDIX  B  -  DEMONSTRATION  PROGRAMS 


/*  This  is  file  "basis.h"  */ 
fdefine     HERMITE      0 
#define     BEZIER       1 
#define     CARDINAL     2 
#define     BSPLINE      3 

/*  the  HERMITE  BASIS  MATRIX  */ 
Matrix  hermitematrix  =  { 

{    2.0,  -2.0,  1.0,  1.0    }, 

{  -3.0,  3.0,-2.0.-1.0  }, 

{0.0,0.0.1.0.0.0      }, 

{  1.0,  0.0.  0.0,  0.0      } 
}: 

/*  the  CARDINAL  BASIS  MATRLX  */ 
Matrix  cardinalmatrix  =  { 

{  -0.5,  1.5,-1.5.  0.5    }, 

{  1.0  ,-2.5,  2.0,-0.5  }, 

{-0.5,0.0,0.5,0.0      }, 

{  0.0,  1.0,  0.0,  0.0       } 

}; 

/*  the  BEZIER  BASIS  MATRLX  */ 
Matrix  beziermatrix  =  { 

{  -1.0,  3.0,-3.0.  1.0    }, 

{  3.0,-6.0,  3.0,  0.0      }, 

{-3.0,3.0,0.0.0.0      }, 

{  1.0,  0.0,  0.0,  0.0       } 

}; 

/*  the  B-SPLINE  BASIS  MATRLX  */ 

Matrix  bsplinematrix  =  { 

{  -1.0/6.0,  3.0/6.0,  -3.0/6.0,  1.0/6.0    }, 
{  3.0/6.0,  -6.0/6.0,  3.0/6.0,  0.0      }, 
{-3.0/6.0,       0.0,3.0/6.0,0.0         }, 
{  1.0/6.0,   4.0/6.0,  1.0/6.0,  0.0      } 

}; 
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/*  This  is  file  "geom.h"  */ 

/*  set  up  the  geometry  matrix  of  x  coordinates  */ 
Coord  geomx[4][4]  =  { 

{  0.0.       100.0,  200.0,  300.0  }, 

{  0.0.       100.0.  200.0.  300.0  }. 

{  1000.0.  900.0,  800.0.  700.0  }, 

{  1000.0.  900.0,  800.0,  700.0  } 
}: 

/*  set  up  the  geometry  matrix  of  y  coordinates  */ 
Coord  geomy[4][4]  =  { 

{  400.0.  500.0.  600.0,  700.0  }, 

{  0.0.     200.0.  400.0,  600.0  }, 

{  0.0.     200.0.  400.0,  600.0  }, 

{  400.0.  500.0,  600.0,  700.0  } 

>; 

/*  set  up  the  geometry  matrix  of  z  coordinates  */ 
Coord  geomz[4][4]  =  { 

{  0.0.  200.0.  400.0,  800.0  }, 

{  0.0.  200.0.  400.0.  800.0  }, 

{  0.0.  200.0,  400.0,  800.0  }, 

{  0.0.  200.0,  400.0,  800.0  } 

>; 
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Program  #1 

This  program  displays  two  wireframe  images  of  the  same 
surface  patch.    The  patch  drawn  in  the  color  YELLOW  is 
produced  by  the  standard  IRIS  patch  functions  and  the 
patch  drawn  in  the  color  RED  is  produced  by  the  parallel 
functions  we  have  developed. 

One  notices  that  the  patch  drawn  via  the  parallel  functions 
has  a  triangular  mesh  appearance  and  that  the  call  to  the 
npatchprecision  does  not  affect  the  displayed  image  as 
does  the  standard  IRIS  patchprecision  function.. 

#include  "gl.h"       /*  IRIS  graphics  library  */ 

^include  "basis. h" 
#include  "geom.h" 

#define     S    CURVES    10 
#define     T    CURVES    10 

main() 
{ 


/*  Loop  variables  */ 
int  pi,  p2; 

/*  initialize  the  graphics  system  */ 

ginit(); 

doublebufferQ: 

gconfig(); 

cursoffQ; 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0.  1023,  0.  767); 

/*  clear  the  graphics  screen  to  BLACK  */ 

color(BLACK); 

clear(); 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
defbasis(BEZIER,  beziermatrix); 
defbasis (CARDINAL,  cardinalmatrix); 
defbasis(BSPLINE,  bsplinematrix) ; 

ndefbasis(BEZIER.  beziermatrix); 
ndefbasis(CARDINAL,  cardinalmatrix); 
ndefbasis(BSPLINE,  bsplinematrix); 
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/**  Specify  how  many  curves  in  each 

parametric  direction.  **/ 
patchcurves(S   CURVES  .  T   CURVES): 
npatchcurves(S   CURVES  .  T   CURVES); 

/**  Make  the  basis  matrices  different 

for  each  parametric  direction.  **/ 
patchbasis(BEZIER.  CARDINAL); 
npatchbasis(BEZIER,  CARDINAL); 

/**  Cycle  through  the  patch  changing  the 

precision  that  the  individual  curves 

comprising  the  patch  are  drawn.  **/ 
for(pl  =  10,  p2  =  100;  pi  <  100;  pi  +  =  5,  p2  -=  5)  { 

/*  Draw  the  image  via  the  IRIS  functions  */ 

viewport(0.  511,  0.  767); 

color(BLACK); 

clear(): 

color(YELLOW); 

patchprecision(pl,  p2): 

patch (geomx,  geomy,  geomz); 

/*  Draw  the  image  via  the  parallel  functions  */ 

viewport(512,  1023,  0.  767); 

color(BLACK); 

clear(): 

color  (RED); 

npatchprecision(pl.  p2); 

npatch(geomx,  geomy,  geomz); 

/*  display  the  wireframe  images  */ 

swapbuffers(); 

sleep(l); 

/*  clear  the  screen  */ 

color(BLACK); 

clear(); 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color(BLACK); 

clear(); 

gexit(); 
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Program  #2 

This  program  illustrates  how  the  triangles  formed  during 
a  surface  patch  decomposition  can  be  put  into  an  IRIS 
graphical  object  and  subsequently  displayed. 

#include  "gl.h"        /*  IRIS  graphics  library  */ 

^include  "basis. h" 
jj include  "geom.h" 

#define     ON  1 

#define     S    CURVES    10 
#define     T   CURVES    10 

/*  Where  we  put  our  intercepted  triangles  */ 
Object  intercepted   object; 

main() 
{ 


/*  declare  the  intercept  function  */ 
int  intercept   functionQ; 

/*  initialize  the  graphics  system  */ 

ginit(); 

doublebuffer  ( ) ; 

gconfig(); 

cursoff(); 

/*  Make  the  intital  object.  */ 

makeobj (intercepted   object  =  genobj()); 

closeobj(); 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport (0,  1023,  0,  767); 

/*  clear  the  graphics  screen  to  BLACK  */ 

color(BLACK); 

clear(); 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis(C ARDINAL,  cardinalmatrix) ; 
ndefbasis(BSPLINE,  bsplinematrix) ; 

/**  Specify  how  many  curves  in  each 

parametric  direction.  **/ 
npatchcurves(S   CURVES  ,  T   CURVES); 
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/**  Make  the  basis  matrices  different 

for  each  parametric  direction.**/ 
npatchbasis(BEZIER.  BSPLINE); 

/**  Set  up  an  intercept  function  to  grab 

the  triangles  from  the  surface  patch 

and  turn  it  on.  **/ 
Set    UserRoutine   for   npatch  (intercept   function); 
User   Routine(ON)~: 

/*  Call  the  npatch  function.  */ 
npatch(geomx.  geomy,  geomz); 

/*  Display  the  surface  patch.  */ 

color(YELLOW); 

callobj (intercepted   object); 

swapbuffers  ( ) ; 

sleep(lO); 

/*  clear  the  graphics  screen  and  exit  */ 

color(BLACK): 

clear(); 

gexit(); 


-a***************************************************** 

This  is  the  user-defined  intercept  function  that 
handles  the  individual  triangles  generated 

during  the  decomposition  of  a  surface  patch. 

♦•♦♦a**************************************************/ 

int  intercept   function  (triangle) 
float  triangle [S] [3]; 


{ 


/*  Open  up  the  object  and  put  in  the  triangles  */ 
editobj(interceptedobject); 
move(triangle[0][0],~triangle[0][l].  triangle [0] [2]); 

draw(triangle[l][0],  triangle[l][l],  triangle[l][2]) 
draw(triangle[2][0],  triangle[2][l],  triangle[2][2]) 
draw(triangle[0][0],  triangle[0][l],  triangle[0][2]) 
closeobj(); 
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Program  #3 

This  programs  illustrates  how  the  intercept  function 
can  be  changed  as  the  program  runs. 

a************************************************************/ 

#include  "gl.h"        /*  IRIS  graphics  library  */ 

finclude  "basis. h" 
#include  "geom.h" 

/**  Define  a  new  type  that  is  a  pointer  to  a  function 

that  returns  an  integer.  **/ 
typedef    int  (*Function   pointer) ()  ; 

#define     ON  1 

fdefine     S   CURVES    10 
#define     T   CURVES    10 

main() 

{ 


/*  Loop  variable  */ 
int  count; 

/*  Declare  an  array  of  pointers  to  functions  */ 
Function   pointer  intercept    functions[4]; 

/*  declare  the  intercept  functions  */ 
int  intercept   functionl(): 
int  intercept_function2(); 
int  intercept    functions ( ) : 
int  intercept_function4(); 

/**  Initialize  the  array  of  pointers  to  the 

intercept  functions.  **/ 
intercept   functions[0]  =  intercept   functionl; 
intercept   functions[l]  =  intercept  function2; 
intercept_functions[2]  =  intercept    functions: 
intercept   functions [3]  =  intercept   function-!: 

/*  initialize  the  graphics  system  */ 

ginit(); 

doublebuffer  ( ) ; 

gconfig(); 

cursoffQ; 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 
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/*  clear  the  graphics  screen  to  BLACK  */ 

color(BLACK); 

clear(); 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis( CARDINAL,  cardinalmatrix); 
ndefbasis(BSPLINE.  bsplinematrix); 

/**  Specify  how  many  curves  in  each 

parametric  direction.  **/ 
npatchcurves(S_CURVES  .  T   CURVES); 

/**  Make  the  basis  matrices  different 
for  each  parametric  direction.  **/ 
npatchbasis(BSPLINE,  CARDINAL); 

/*  Initially  use  the  default  intercept  function  */ 
Set    Default    Routine   for   npatch(); 
User   Routine  (ON): 

/*  Step  through  each  intercept  function  */ 
for(count  =  0;  count  <  4  ;  countH — h)  { 

/**  Set  up  an  intercept  function  to  grab 

the  triangles  from  the  surface  patch.  **/ 
Set    User   Routine   for  npatch(intercept_functions[count]); 

/**  Call  the  npatch  function  using  the  current 

intercept  function.  **/ 
npatch(geomx,  geomy,  geomz); 

/*  Display  what  the  intercept  function  did.  */ 

swapbuffers  ( ) ; 

sleep(2); 

/*  Clear  the  screen  and  do  another  one.  */ 

color(BLACK); 

clear(); 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color  (BLACK); 

clear(); 

gexit(); 
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,'*******************************************************» 

This  intercept  function  draws  each  triangle  RED. 

******************************************************** 

int  intercept    function  1  (triangle) 
float  triangle]3][3]; 

{ 

color(RED); 
poly (3,  triangle): 

} 


/******************************************************** 

This  intercept  function  draws  each  triangle  YELLOW. 

******************************************************** ^ 

int  intercept   function2 (triangle) 
float  triangle[3][3]; 

{ 

color(YELLOW): 

poly (3,  triangle): 
} 


/******************************************************** 

This  intercept  function  draws  each  triangle  GREEN. 

**********************••********************••**********/ 

int  intercept   function3 (triangle) 
float  triangle[3][3]; 

{ 

color  (GREEN): 

poly(3,  triangle): 
} 


/*****************************••************************* 

This  intercept  function  draws  each  triangle  BLUE. 

***************************************•****************/ 

int  intercept   function4 (triangle) 
float  triangle[3][3]; 


color  (BLUE); 
poly (3,  triangle); 
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Program  #4 

This  program  illustrates  how  the  intercept  function 
can  be  combined  with  an  illumination  model  to 
provide  a  realistic,  illuminated,  three-dimensional 
solid-filled  curved  surface. 

•♦-A********************************************************/ 

#include  "gl.h"       /*  IRIS  graphics  library  */ 

^include  "device. h" 
^include  "basis. h" 
^include  "geom.h" 

#define    ON  1 

#define    S    CURVES    25 
#define    T   CURVES    25 

main() 
{ 


/*  Declare  the  intercept  function.  */ 
int    light   polyQ; 

/*  loop  variables  */ 
int   ij; 

/*  Initialize  the  graphics  system.  */ 

ginitQ; 

singlebuffer(); 

gconfig(); 

cursoffQ; 

/*  Clear  the  display.  */ 

color(BLACK); 

clear(); 

/*  Set  up  new  viewing  parameters.  */ 
ortho(0.0,  1023.0.  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0.  1023,  0,  767); 

/*  Clear  drawing  area.  */ 

color(CYAN); 

clear(); 

/*  Set  for  hidden  surface  elimination.  */ 

setdepth(0x3FFF,  OxCOOO); 

zclear(); 

zbuffer(TRUE); 
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/*  Load  the  color  map  ramp  with  a  grey  scale.  */ 
for  (i  =  0;i  <    256;  i++) 

{ 

mapcolor(8-fi.  i,  i,  i); 

} 

/*  Associate  an  id  with  a  basis  matrix.  */ 
ndefbasis(BEZIER,  beziermatrix): 
ndefbasis(CARDINAL,  cardinalmatrix); 
ndefbasis(BSPLINE,  bsplinematrix); 

/**  Provide  a  different  basis  for  each 

parametric  direction.  **/ 
npatchbasis(CARDINAL.  BEZIER); 

/**  Provide  the  number  of  curves  in 

each  parametric  direction.  **/ 
npatchcurves(S   CURVES  .  T   CURVES); 

/*  Set  up  the  intercept  function.  */ 

Set    UserRoutinefornpatch  (light   poly): 

U ser   Routine  ( O N ) ; 

while(TRUE) 

{ 

/*  Clear  the  z  buffer.  */ 
zclear(); 

/*  Hold  display  if  MOUSE2  is  down.  */ 
if(getbutton(MOUSE2)) 

{ 

/*  Resume  when  MOUSEl  is  pressed  */ 

while(!getbutton(MOUSEl))  ; 
} 

/*  Exit  when  MOUSEl  and  MOUSE2  and  MOUSE3  are  down.  */ 
if(getbutton(MOUSEl)  &&  getbutton(MOUSE2)  &&  getbutton(MOUSE3)) 
break; 

/*  Clear  the  drawing  area.  */ 

color  (CYAN); 

clear(); 

/*  Set  the  current  color.  */ 
color  (BLACK); 

/*  Draw  initial  surface  patch.  */ 
npatch(geomx,  geomy,  geomz); 
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/**  Change  the  y  coordinates  to  get  a 

different  surface  patch  the  next 

time  we  draw  display.  **/ 
for(i=l;  i<3;i++) 

{ 

for(j=0;j<4;j+  +  ) 

{ 

geomy[i][j]  =  (float) (lrand48()  %  900); 

} 
} 

/*  Draw  the  surface  patch.  */ 
npatch(geomx,  geomy,  geomz): 

/**  Change  the  y  coordinate  values  to 

make  another  surface  patch.  **/ 
for(i=l;  i<3;  i++) 

{ 

for(j=0:  j<4:j+  +  ) 

{ 

geomy[i][j]  =  (lrand48()%2)*geomy[i]  [j] 

} 
} 
} 

/*  Clean  up  and  exit  the  program.  */ 

color(BLACK); 

clear(): 

gexit(); 


The  user-defined  intercept  function  used  to  grab  the 
triangles  generated  in  the  surface  patch  decomposition. 

int  light   poly(Triangle) 
Coord  Triangle[3][3]; 

{ 

/*  Put  each  triangle  through  an  illumination  model.  */ 

lightthepoly (Triangle,  3,  350.0,  -1750.0,  350.0,  350.0,  1750.0,  350.0,  9,  264); 


} 
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lightpoly.c 

It  is  a  routine  that  computes  lighting  for  a  polygon  based 
upon  the  angle  between  the  Normal  vector  of  the  polygon 
and  the  direction  to  the  light  source. 

lightthepoly(xyz,ncoords,ax,ay,az,lx,ly,lz,colormin,colormax) 

xyz[][3]  =  floating  coords  of  the  polygon. 

ncoords  =  number  of  coordinates. 

ax.ay.az  =  interior  point  of  the  ■whole  object.  Used  to  determine 
outward  facing  normal  of  the  polygon.  This  is  the  same 
point  of  reference  that  would  be  used  for  backface 
polygon  removal. 

lx,ly.lz  =  vector  pointing  in  direction  of  the  light  source. 

colormin.  colormax  =  indices  used  for  the  colors  assigned  to  this 
polygon.  The  user  is  responsible  for  setting 
up  the  color  ramp. 

Note:  the  routine  also  puts  the  polygons  out  ordered  counterclockwise 
with  respect  to  the  interior  point  for  ease  of  backface  polygon 
removal. 

#include  <gl.h> 
#include  <math.h> 

fdefine  PIDIV2        1.570796327 
fdefine  CLOCKWISE   1 
#define  ROW  3 

lightthepoly(xyz,ncoords,ax,ay,az,lx,ly,lz,colormin.colormax) 

Coord  xyz[][3]; 

unsigned  int   ncoords; 

Coord  ax.ay,az:     /*  interior  point  of  the  whole  object.  */ 

Coord  lx,ly,lz;     /*  direction  to  the  light  source  */ 

int  colormin,colormax;     /*  color  min  /max  indices         */ 

{ 
/*  temp  coord  hold  */ 
Coord  *txyz; 

/*  loop  temps  */ 

register  unsigned  short  int     i  j; 
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/*  direction  test  function  */ 
int  npoly   orient  (); 

/*  vectors  used  to  compute  the  polygon's  normal  */ 
Coord  vl[3].v2[3]: 

/*  the  polygon's  normal  */ 
Coord  normal[3]: 

/*  normal's  magnitude     */ 
Coord  normalmag; 

/*  light's  magnitude     */ 
Coord  lightmag; 

/*  dot  product  of  N  and  L  */ 
double   dotprod: 

/*  angle  between  N  and  L  */ 
float  radians; 

/*  color  to  use  in  drawing  the  polygon  */ 
unsigned  short  int    colortouse: 

/*  allocate  memory  for  a  temporary  array  */ 

txyz  =  (Coord  *)  calloc((ncoords  *  3),  sizeof(Coord)): 

/**  orient  the  polygon  so  that  its  counterclockwise  with  respect 

to  the  interior  point  **/ 
if(npoly  orient  (ncoordspcyz,ax,ay,az)  ==  CLOCKWISE) 

{ 

/*  the  polygon  is  clockwise,  reverse  it.  */ 

for(i=0:  i  <  ncoords;  i=i+l) 

{ 
for(j=0;  j  <  ROW;  j=j+l) 

{ 
*(txyz  +  (i  *  ROW)  +  j)  =  xyz[ncoords-i-l][j]; 

} 

} 

} 
else 

{ 

/*  no  need  to  reverse  */ 
for(i=0;  i  <  ncoords;  i=i+l) 

{ 
for(j=0;  j  <  ROW;  j=j+l) 

{ 
*(txyz+  (ROW*i)+j)=xyz[i][j]; 

} 

} 
} 
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/*  the  coordinates  are  ordered  counterclockwise  in  array  txyz  */ 

/**  compute  the  normal  vector  for  the  polygon  using  the  first 
three  vertices...  **/ 

/*  compute  the  first  vector  to  use  in  the  computation  */ 
vl[0]  =  *(txyz  +  6)  -  *(txyz  +  3);  /*  txyz[2][0]  -  txyz[l][0]  */ 
vl[l]  =  *(txyz  +  7)  -  *(txyz  +  4);  /*  txyz[2][l]  -  txyz[l][l]  */ 
vl[2]  =  *(txyz  +  8)  -  *(txyz  +  5);  /•  txyz[2][2]  -  txyz[l][2]  */ 

/*  compute  the  second  vector  to  use  in  computing  the  normal  */ 
v2[0]  =  *(txyz      )  -  *(txyz  +  3);  /*  txyz[0][0]  -  txyz[l][0]  */ 
v2[l]  =  *(txyz  +  1)  -  *(txyz  +  4);  /*  txyz[0][l]  -  txyz[l][l]  */ 
v2[2]  =  *(txyz  +  2)  -  *(txyz  +  5);  /*  txyz[0][2]  -  txyz[l][2]  */ 

/*  the  normal  is  vl  x  v2  */ 
normal[0]  =  vl[l]*v2[2]  -  vl[2]*v2[l]; 
normaljl]  =  vl[2]*v2[0]  -  vl[0]*v2[2]; 
normal[2]  =  vl[0]*v2[l]  -  vl[l]*v2[0]: 

/*  compute  the  magnitude  of  the  normal  */ 
normalmag  =  sqrt((normal[0]*normal[0])  + 

(normal  [l]*normal[l] ) + 

(normal[2]*normal[2])); 


/*  compute  the  magnitude  of  the  light  */ 
lightmag  =  sqrt((lx  *  lx)  +  (ly  *  ly)  +  (lz  * 


i«)); 


/*  compute  N  .  L  (normal  dot  product  with  the  light  source  direction)  */ 
dotprod  =  (double)  ((normalfO]  *  lx) 

+  (normaljl]  *  ly) 

+  (normal[2]  *  lz)); 

/*  compute  the  unit  normal  */ 

dotprod  =  (double)  ((dotprod/ (normalmag  *  lightmag))); 

/*  dotprod  =  cos(theta)  of  the  angle  between  N  and  L. 

Convert  to  angle  in  radians  */ 
radians  =  acos(dotprod); 
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/*  compute  the  color  we  should  use  */ 
if(-PIDIV2  <=  radians  &&  radians  <=  PIDIV2) 

{ 
/*  if  the  angle  is  negative,  set  to  positive  */ 
if(radians  <  0.0) 

{ 
radians  =  -radians; 

} 

colortouse  =  ((colormax-colorirun)/PIDIV2)*(PIDIV2-radians)+colormin: 

} 
else 

{ 
colortouse  =  colormin; 

} 

/*  set  the  color  */ 
color  (colortouse); 

/*  draw  the  poly  */ 
polf(ncoords.txyz) : 

/*  free  up  memory  allocated  for  the  temporary  array  */ 
cfree(txyz.  (ncoords  *  3).  sizeof( Coord)); 


} 


55 


:*****♦*************************************    +    +     +     *     +    +    *********** 

npoly  orient. c 

This  routine  determines  a  polygon's  orientation 
with  respect  to  its  normal  and  a  reference  point. 
Orientation  is  either  clockwise  or  counter-clockwise. 
The  point  of  reference  must  not  lie  in  the  polygon's 
plane. 

#include  <gl.h> 
^include  <math.h> 

int  npoly  orient  (ncoords.xyz.xinside,yinside,zinside) 

unsigned  int   ncoords; 

Coord  xyz[][3]; 

Coord  xinside.  yinside.  zinside: 

{ 
/*  loop  temps  */ 

register  unsigned  short  int     i  j: 

/*  center  coordinate  of  the  polygon  */ 
Coord  center [3]; 

/**  vector  hold  locations  for  the  vectors  that  run 
from  the  center  coordinate  to  the  points  of  the 
polygon  **/ 

Coord  a[3],b[3]; 

/**  points  on  line  containing  normal  that  are 

on  opposite  sides  of  the  plane  containing 

the  polygon.  **/ 
Coord  xn[3],  xmn[3]; 

/*  distance  to  point  n  from  pt  inside.  */ 
float  distton; 

/*  distance  to  point  -n  from  pt  inside.  */ 
float  disttomn; 

/*  the  normal  vector  computed  from  a  x  b  */ 
Coord  normal  [3]; 

/*  compute  the  center  coordinate  of  the  polygon  */ 
center[0]  =  0.0; 
centerjl]  =  0.0 
center  2    =  0.0 
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for(i=0;  i  <  ncoords;  i-\ — (-) 

{ 

for(j=0;j  <  3;j+  +  ) 

{ 
center[j]  +  =  xyz[i][j]; 

} 
} 

/*  divide  out  by  the  number  of  coordinates  */ 
for(j=0;j  <  3;j++) 

{ 
center[j]  =  center[j]/(float)ncoords; 

} 

/**  check  the  first  2  coordinates  of  the 
polygon  for  their  direction  **/ 

/**  compute  vector  a.  It  runs  from  the 
center  coordinate  to  coordinate  0  **/ 
for(j=0;j  <  3:j+  +  ) 

{ 
a[j]  =  xyz[0][j]  -  center[j]; 

} 

/**  compute  vector  b.  It  runs  from  the 
center  coordinate  to  coordinate  1  **/ 
for(j=0;j  <3;j+  +  ) 

{ 

b[j]  =  xyz[l][j]  -  center[j]; 

} 

/*  compute  a  x  b  to  get  the  normal  vector  */ 
normal[0]  =  a[l]*b[2]  -  a[2]*b[l] 
normaljl]  =  a[2]*b[0]  -  a[0]*b[2] 
normal[2]  =  a[0]*b[l]  -  a[l]*b[0] 

/**  compute  point  n,  offset  pt  from  center  in 

direction  of  normal  **/ 
for(j=0:j<3;j++) 

{ 
xn[j]  =  center[j]  +  normal[j]; 

} 

/**  compute  point  -n,  offset  pt  from  center 

in  opposite  direction  from  normal.  **/ 
for(j=0;j  <  3;j++) 

{ 
xmn[j]  =  center  [j]  -  normal [j]; 

} 
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/*  compute  the  distance  the  inside  pt  is  from  point  n  */ 
distton  =  sqrt((xn[0]  -  xinside)  *  (xn[0]  -  xinside)  -f 

(xn[l]  -  yinside)  *  (xn[l]  -  yinside)  + 

(xn[2]  -  zinside)  *  (xn[2]  -  zinside)); 

/*  compute  the  distance  the  inside  pt  is  from  point  -n  */ 
disttomn  =  sqrt((xmn[0]  -  xinside)  *  (xmn[0]  -  xinside)  + 
(xmn[l]  -  yinside)  *  (xmn[l]  -  yinside)  ■+ 
(xmn[2j  -  zinside)  *  (xmn[2]  -  zinside)); 

/**  if  the  dist(n)  <  dist(-n).  then  n  points  back  towards  the 
inside  point  and  is  on  the  same  side  of  the  plane  as  inside, 
a  x  b  is  then  clockwise.  **/ 

if(distton  <  disttomn) 

{ 
return(l);     /*  clockwise  */ 

} 
else 

{ 
return(O);     /*  counterclockwise  */ 

} 


} 


58 


APPENDIX  C         BENCHMARK  PROGRAM 


********■******  +  ************************************  +  **♦****** 

This  is  the  BENCHMARK  PROGRAM  used  to  test  the 
relative  performance  of  the  standard  IRIS  functions 
and  the  parallel  routines  to  those  standard 
functions. 

#include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  "basis. h" 
^include  "geom.h" 

#define   IRIS  /*  Which  set  to  test  switch.  */ 

#define    MAX   TIMES   THRU    100 
#define    S    CURVES  25 

#define    T   CURVES  25 

#define   ONE  1 


main() 

{ 


int  times   thru; 

/*  initialize  the  graphics  system  */ 

ginit(): 

doublebuflFer(): 

gconfig(); 

cursoff(); 

/*  clear  the  graphics  screen  */ 

color(BLACK); 

clear(); 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 

/*  clear  the  graphics  screen  to  CYAN  */ 

color(CYAN); 

clear(); 
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#ifdef  IRIS 

/*  Use  the  standard  IRIS  functions  */ 
defbasis(BEZIER,  beziermatrix); 
defbasis(CARDINAL.  cardinalmatrix); 
defbasis(BSPLINE,  bsplinematrix) : 

patchbasis(BEZIER,  BEZIER); 
patchcurves(S   CURVES  .  T  CURVES); 
patchprecision(ONE,  ONE); 
#else 

/*  Use  the  parallel  functions  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis(C ARDINAL.  cardinalmatrix) ; 
ndefbasis(BSPLINE.  bsplinematrix) ; 

npatchbasis(BEZIER.  BEZIER); 
npatchcurves(S   CURVES  ,  T   CURVES); 
npatchprecision(ONE,  ONE): 
fendif 
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for  (times   thru  =  0:  times   thru  <  MAX   TIMES   THRU:  times   thru++)  { 

/*  clear  the  graphics  screen  to  CYAN  each  time  thru  */ 

color(CYAN); 

clear(); 

/*  draw  the  wireframe  surface  patch  in  BLACK  */ 

color(BLACK): 
#ifdef  IRIS 

/*  Using  the  IRIS  function  */ 

patch(geomx,  geomy.  geomz); 
#else 

/*  Using  the  parallel  function  */ 

npatch(geomx,  geomy,  geomz); 
#endif 

/*  display  the  wireframe  image  */ 

swapbuffers  ( ) ; 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color(BLACK); 

clear(); 

gexit(); 
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APPENDIX  D  -  PARALLEL  FUNCTIONS  SOURCE  CODE 


/a*********************************************************** 

FILE  npatch.c 

Author  Gary  W.  TAYLOR  (Captain  USMC) 

Date  1  December  1986 

Place  Naval  Postgraduate  School,  Monterey  CA 

Environment  


Silicon  Graphics.  Inc.,  IRIS  2400 
graphics  workstation,  UNIX  operating 
system  (GL2-W3.4). 


Purpose 


The    following  C    source  code   provides  a  set 
of  "shadow"  routines  to  parallel  the  standard 
IRIS   2400  graphics  workstation  surface  patch 
routines.     These   parallel  routines   provide 
their   user   the   capability   to   generate 
solid-filled  parametric  bicubic  surface  patches. 


Notes 


As  of  the  current  date,  there  are  no  known 
side-effects  or  bugs  associated  with  using 
these  functions. 


Limitations  

These  functions  can  be  used  in  immediate  mode  only. 

ft*********************************************************** 
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#include  "gl.h"       /*  IRIS  graphics  library  */ 
#include  "stdio.h" 

#define    OFF  0 

/*  Here  is  how  we  define  a  surface  point.  */ 
typedef  struct  { 

Coord  x: 

Coord  y; 

Coord  z: 
}  Point: 

/*  Structure  used  to  track  user  supplied  basis  matrices.  */ 
static  struct  list    elem  { 

/*  Matrix  id  number.  */ 

long   nid: 

/*  Pointer  to  the  basis  matrix.        */ 
float    *nmatrix; 

/*  Pointer  to  the  next  basis  matrix.  */ 
struct  list   elem  *next_elem_ptr; 
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/*  Used  in  forward  difference  computation  in  the  U  direction.  */ 
static    Matrix  Precision   Matrix  U; 

/*  Used  in  forward  difference  computation  in  the  V  direction.  */ 
static    Matrix  PrecisionMatrixV; 

/*  The  default  intercept  function.  */ 
static  int      Default   systemroutine  (); 

/*  Pointer  to  the  currently  defined  U  basis  matrix.  */ 
static  float  *current    Ubasis; 

/*  Pointer  to  the  currently  defined  V  basis  matrix.  */ 
static  float  ""current    Vbasis; 

/*  How  many  curves  in  the  U  direction.  */ 
static  long  UCURVES: 

/*  How  many  curves  in  the  V  direction.  */ 
static  long  VCURVES; 

/*  How  may  curve  segments  in  the  U  direction.  */ 
static  long  USEGMENTS; 

/*  How  may  curve  segments  in  the  V  direction.  */ 
static  long  VSEGMENTS; 

/*  Pointer  to  linked  list  of  user  supplied  basis  matrices.*/ 
static  struct  listelem  *head_of  list  =  NULL; 

/*  Set  initial  user  routine  to  the  default.  */ 

static  int    (*   User  routine)  ()  =     Default   system  routine; 

/*  Initially  do  not  call  user's  function.  */ 
static  int      User   routine  is   on  =  OFF; 
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Default    system   routine 

This  routine  is  -what  the  system  automatically  does  if  the 
user  does  not  supply  a  particular  intercept  function. 

******  41414141  ************************  414141  4c*  ************************* 

static  int      Default   system  routine  (Triangle) 

Coord  Triangle [3]  [3]; 

{ 

/*  Draw  the  polygon  with  a  standard  IRIS  function.  */ 
poly  (3,  Triangle); 

} 


/**************************************************************** 

Set    User   Routinefornpatch 

This  function  allows  the  user  to  supply  an  intercept  function 
to  be  used  to  manipulate  the  triangles  that  are  generated  in 
the  decomposition  of  a  surface  patch. 

****************************************************************/ 

void  Set   User   Routine  for  npatch  (routine) 

int     (""routine)  (); 

{ 

/*  Save  the  pointer  to  the  user's  function.  */ 
User   routine  =  routine; 

} 
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User   Routine 

This  function  allows  the  user  to  turn  the  intercept  function 
on  or  off  at  will. 

*********+****************+*+***************♦*****+*++*********+ 

void  User   Routine  (boolean) 
int     boolean; 

{ 

If  boolean    =  0  then  the  intercept  function  is  turned  OFF. 
If  boolean  !=  0  then  the  intercept  function  is  turned  ON. 

V 

User  routine   is   on  =  boolean: 


/it*************************************************************** 

Set    Default   Routinefor   npatch 

This  function  allows  the  user  to  reset  the  intercept  routine 
to  the  same  routine  used  by  the  system. 

void  Set   Default   Routine  for  npatch  () 

{ 

/*  Point  to  the  default  intercept  function.  */ 
User   routine  =     Default   system  routine; 


} 
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ndefbasis 

This  function  is  equivalent  to  the  IRIS  defbasis  function 
in  that  it  allows  the  user  to  define  a  basis  matrix  for 
use  in  the  generation  of  patches,   matrix  is  saved  and 
is  associated  with  id.    id  may  then  be  used  in  subsequent 
calls  to  npatchbasis. 

void  ndefbasis  (id,  matrix) 
long   id: 
Matrix  matrix; 

{ 

/*  Special  processing  first  time  this  function  is  called.  */ 
static  int    has   been   called  =  FALSE; 

/*  Data  structure  pointer  for  new  entry.  */ 
struct  list    elem  *new   elem  to  add; 

/*  Pointer  to  search  linked  list.  */ 
struct  list    elem  *walking_ptr; 

/*  Pointer  to  a  copy  of  the  user  supplied  basis  matrix.  */ 
float    *pmatrix: 

/*  Loop  variables  */ 
int     row: 
int     column; 

/*  Get  memory  for  the  new  data  elements.  */ 

new   elem   to   add  =  (struct  list   elem  *)  malloc  (sizeof  (struct  list_elem)); 
pmatrix  =  (float  *)  calloc  (sizeof  (Matrix),  sizeof  (float)); 

/*  Make  a  copy  of  the  basis  matrix  passed  in  by  the  user.  */ 
for  (row  =  0:  row  <  4;  rowH — \-) 
for  (column  =  0:  column  <  4;  column++)  { 

*  (pmatrix  +  (4  *  row)  +  column)  =  ma trixfrow]  [column]; 

} 

/*  Associate  the  user  supplied  id  to  this  basis  matrix.  */ 
new   elem   to   add  ->  nid  =  id; 

/*  Point  to  the  copied  basis  matrix.  */ 
new   elem  to   add  ->  nmatrix  =  pmatrix; 
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/**  Determine  how  to  add  this  information  into  the  linked  list 
of  basis  matrices.  **/ 

/*  Does  a  list  already  exist?  */ 
switch  (has   been   called)  { 

case  TRUE: 

/*  Point  to  the  beginning  of  the  list.  */ 
walking   ptr  =  head   of  list; 

/**  Traverse  the  list  looking  to  see  if  the  id  number 
already  exists.  **/ 

while  ((walking   ptr  ->  nid  !=  id)  && 
(walking   ptr  ->  next   elem  ptr  !=  headoflist))  { 

/*  Walk  through  the  linked  list.  */ 

walking   ptr  =  walking  ptr  ->  next   elem  ptr; 

} 

/*  id  already  exists  so  we  can  reuse  its  allocated  memory.  */ 
if  (walkingptr  ->  nid  ==  id)  { 

/*  Get  rid  of  the  old  basis  matrix.  */ 

cfree  (walking   ptr  ->  nmatrix,  sizeof  (Matrix),  sizeof  (float)); 

/*  Point  to  the  replacement  matrix.  */ 
walking  ptr  ->  nmatrix  =  pmatrix; 

/*  Get  rid  of  the  un-needed  data  structure.  */ 

cfree  (new   elem   to   add,  1,  sizeof  (struct  list   elem)); 

} 

else  {  /*  The  id  does  not  exist    */ 

/**  Manipulate  the  pointers  to  add  the  new 

data  element  to  the  linked  list.  **/ 
newelem  to  add  ->  next   elem  ptr  =  head   of  list  ->  next   elem  ptr: 
headoflist  ->  next   elem  ptr  =  new   elem   to   add; 

} 
break; 


68 


case  FALSE: 

/*  No  linked  list  of  basis  matrices  exists  so  we  start  up  one.  */ 

/*  Create  the  pointer  to  the  front  of  the  list.  */ 
head_of_list  =  new  elem  to  add; 
head   of  list  ->  next   elem   ptr  =  head   of  list: 

/*  Make  sure  we  do  not  do  this  again.  */ 
has   been   called  =  TRUE: 

break: 

default: 
break; 

} 


} 


69 


npatchbasis 

This  function  is  equivalent  to  the  IRIS  patchbasis  function 
in  that  it  sets  the  current  basis  matrices  for  both  the 
U  and  V  parametric  directions  of  a  surface  patch.    The  current 
U  and  V  bases  are  used  when  the  npatch  command  is  issued. 

***************+*****#********+*******+***+*»*+***♦+************ 

int     npatchbasis  (uid,  vid) 

long   uid,  vid; 

{ 

struct  list   elem  *walking_ptr; 

/*  ERROR:  no  linked  list  of  basis  matrices.  */ 
if  (head   of  list  ==  NULL)  { 

fprintf  (stderr,  "Opatchbasis:  no  basis  matrices  definedO); 
exit  (-1); 

} 

/*  Traverse  the  list  looking  for  the  desired  U  basis   matrix.  */ 
walking   ptr  =  head   of  list; 
while  ((walking   ptr  ->  nid  !=  uid)  && 

(walking   ptr  ->  next   elem   ptr  !=  head   of  list))  { 

walking   ptr  =  walking   ptr  ->  next   elem  ptr; 

} 

if  ((walkingptr  ->  nid  !=  uid)  && 

(walkingptr  ->  next   elemptr  ==  head   of  list))  { 

/*  ERROR:  U  basis  matrix  does  not  exist  in  the  linked  list.  */ 

fprintf  (stderr,  "Opatchbasis:  undefined  U  basis  matrix  %d0,  uid); 
exit  (-1); 

} 

current   U   basis  =  walking  ptr  ->  nmatrix; 
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/*  Traverse  the  list  looking  for  the  desired  V  basis   matrix.  */ 
walking   ptr  =  head   of  list: 
while  ((walking   ptr  ->  nid  !=  vid)  &&: 

(walking   ptr  ->  next    elem  ptr  !=  head   of  list))  { 

walking   ptr  =  walking   ptr  ->  nextelemptr; 

} 

if  ((walking   ptr  ->  nid  !=  vid)  && 

(walking   ptr  ->  nextelemptr  ==  head   of  list))  { 

/*  ERROR:  V  basis  matrix  does  not  exist  in  the  linked  list.  */ 

fprintf  (stderr,  "Opatchbasis:  undefined  V  basis  matrix  %d0.  vid): 
exit  (-1): 

} 

current    V   basis  =  walking   ptr  ->  nmatrix; 
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/************************#*************+*****+******++****+*****+ 

npatchcurves 

This  function  is  similar  to  the  IRIS  patchcurves  command. 
ucurves  and  vcurves  set  the  subdivision  parameters  used  in 
decomposing  the  surface  patch.    An  individual  patch  will  be 
decomposed  into  a  (ucurves  -1)  *  (vcurves  -1)  grid  with 
each  grid  generating  two  triangular  polygons. 

****************************************************************; 

void  npatchcurves  (ucurves,  vcurves) 

long    ucurves,  vcurves; 

{ 

/*  Prevent  an  inappropriate  number  of  curves.  */ 
UCURVES  =  (ucurves  <  2  ?  2  :  ucurves  -  1); 
VCURVES  =  (vcurves  <  2  ?  2  :  vcurves  -  1); 

/**  Set  up  the  Precision   Matrix   U  used  in  the  forward  difference 
along  the  U  direction.  **/ 


Precision 

Matrix 

U[0] 

[0] 

=  6.0  /  (float) 

Precision 

Matrix 

U[l] 

[o] 

=  6.0  /  (float) 

Precision 

Matrix 

U[l] 

[1] 

=  2.0  /  (float) 

Precision 

Matrix 

U[2] 

[0] 

=  1.0  /  (float) 

Precision 

Matrix 

~U[2] 

[1] 

=  1.0  /  (float) 

Precision 

Matrix 

~U[2] 

[2] 

=  1.0  /  (float) 

Precision 

Matrix 

~U[3] 

[3] 

=  1.0; 

(UCURVES  *  UCURVES  *  UCURVES); 
(UCURVES  *  UCURVES  *  UCURVES); 
(UCURVES  *  UCURVES); 
(UCURVES  *  UCURVES  *  UCURVES); 
(UCURVES  *  UCURVES); 
(UCURVES): 


/**  Set  up  the  Precision   Matrix   V  used  in  the  forward  difference 
along  the  V  direction.    **/ 


Precision 
Precision 
Precision 
Precision 
Precision 
Precision 
Precision 


Matrix 
Matrix 
Matrix 
Matrix 
Matrix 
Matrix 
Matrix 


V[0][0] 
V[1][0] 

v[i][i] 

V[2][0] 
V[2][l] 
V[2][2] 
V[3][3] 


=  6.0  /  (float 
=  6.0  /  (float 
=  2.0  /  (float 
=  1.0  /  (float 
=  1.0  /  (float 
=  1.0  /  (float 
=  1.0; 


(VCURVES  *  VCURVES  *  VCURVES); 
(VCURVES  *  VCURVES  *  VCURVES); 
(VCURVES  *  VCURVES); 
(VCURVES  *  VCURVES  *  VCURVES): 
(VCURVES  *  VCURVES); 
(VCURVES); 


} 
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npatchprecision 

This  function  is  similar  to  the  IRIS  patchprecision  routine 
but  is  only  used  to  maintain  consistency  with  the  IRIS 
routines.    Its  results  are  not  used  by  any  other  function 
in  the  suite. 

it**************************************************************/ 

void  npatchprecision  (usegments.  vsegments) 

long   usegments.  vsegments; 

{ 

/*  Prevent  an  inappropriate  number  of  segments.  */ 
USEGMENTS  =  (usegments  <  2  ?  2  :  usegments  -  1); 
VSEGMENTS  =  (vsegments  <  2  ?  2  :  vsegments  -  1); 

} 
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nspeckle 

This  function  finishes  computation  and  hands  off  each  triangle 
to  an  intercept  function. 

♦♦♦a************************************************************/ 

static    void  nspeckle  (Coord   array) 
Point  *  Coord   array [4]; 

{ 

register    Point  *  Patch   array; 

/*  Get  enough  memory  to  hold  all  the  points  that  will  be  generated.  */ 
Patch   array  =  (Point  *)  calloc  ((UCURVES  +  1)  *  (VCURVES  +  1), 
sizeof  (Point)); 

{ 

/*  Pointer  to  a  particular  point.  */ 
register    Point  *  Where; 

register  unsigned  int  total   points; 
register  unsigned  int  point    count; 
register  unsigned  int  t   count; 
register  unsigned  int  count; 
register  unsigned  int  Row; 
register  unsigned  int  Column; 

/*  Used  in  generating  points  on  the  surface.  */ 
Matrix  control  matrix; 

/*  Intermediate  matrix  to  hold  mathematical  results.  */ 
Matrix  interl; 

/*  For  every  point  in  the  the  U  parametric  direction.  */ 
for  (pointcount  =  0,  total   points  =  0; 

point   count  <=  UCURVES;  point   count++)  { 

/*  Build  a  control  matrix  for  the  current  curve.  */ 
for  (count  =  0;  count  <  4;  count++)  { 

Where  =  (Point  *)  (Coord   array  [count]  +  pointcount); 

control_matrix[count][0]  =  Where  ->  x; 
control_matrix[count][l]  =  Where  ->  y; 
control_matrix[count][2]  =  Where  ->  z; 

} 
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/**  Generate  the  matrix  to  compute  the  forward  difference  on. 
The  forward  difference  matrix  is  equal  to 
PrecisionMatrixV  x  current   V   basis  x  control   matrix.  **/ 

pushmatrix  (); 

loadmatrix  (control   matrix); 

multmatrix  (current    V   basis): 

multmatrix  (Precision   Matrix  V); 

getmatrix  (interl); 

popmatrix  (); 

/*  Generate  the  points  on  the  curve  in  the  V  direction  */ 
for  (t   count  =  0;  t   count  <=  VCURVES:  t   count  +  +  .  total   points+  +  )  { 

(Patch  array  +  total  points)  ->  x  =  interl[3][0j: 
(Patcharray  +  total  points)  ->  y  =  interl[3][l]; 
(Patch   array  +  total   points)  ->  z  =  interl[3][2]; 

/*  Do  the  forward  difference.  */ 
for  (Row  =  3;  Row  >  0:  Row--) 
for  (Column  =  0:  Column  <  4;  Column-f+) 
interl  [Row]  [Column]  =  interl  [Row]  [Column]  +  interl  [Row  -  1]  [Column] 

} 
} 
} 
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{ 


/*  Place  to  put  a  triangle  to  send  out  to  user  */ 
Coord  Triangle    1  [4]  [3]; 

register  Point  *  Where; 
register  unsigned  int  Row; 
register  unsigned  int  Column; 

/*  Decompose  the  patch  into  its  individual  triangles.  */ 
for  (Row  =  0:  Row  <  UCURVES;  Row++) 
for  (Column  =  0;  Column  <  VCURVES;  Column+  +  )  { 

Where  =  (Patch   array  +  ((VCURVES  +  1)  *  Row)  +  Column); 


Triangle  1[0 

Triangle  l[0 

Triangle  l[0 

Triangle  1[1 

Triangle  1[1 

Triangle  l[l 

Triangle  1  [2 

Triangle  1  [2 

Triangle  1  [2 

Triangle  1  [3 

Triangle  1  [3 

Triangle  13 


0]  =  Where  ->  x; 
1]  =  Where  ->  y: 
2]  =  Where  ->  z; 

01  =  (Where  +  1) 
1]  =  (Where  +  1) 

2  =  (Where  +  1 


>  x; 

>  y: 

->  z; 


0]  = 

1]  = 
2}  = 

0]  = 

1]  = 
2]  = 


(Where  +  VCURVES  +  1)  ->  x; 
(Where  +  VCURVES  +  1)  ->  y; 
(Where  +  VCURVES  +  1)  ->  z; 


■>  x; 


(Where  +  VCURVES  +  2) 
(Where  +  VCURVES  +  2)  ->  y; 
(Where  +  VCURVES  +  2)  ->  z; 


} 


/*  Does  the  user  have  an  intercept  routine?  */ 
if  (Userroutineis   on)  { 
/*  Yes*/ 

(*    User   routine)  (fcTriangle _l[0][0]); 
(*   User  routine)  (fcTriangle   l[l][0]); 

} 
else  { 

/*  No  user  routine,  so  we  use  the  default.  */ 

Default   system   routine  (fcTriangle    1[0][0]); 

Default   systemroutine  (&Triangle_l[l][0]); 

} 
} 
} 

/*  Return  the  memory  used  to  the  system.  */ 
cfree  (Patch   array,  (UCURVES  +  1)  *  (VCURVES  +  1),  sizeof  (Point)): 
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npatch 

This  function  rearranges  the  input  matrices  into  a  form  readily 
used  in  computing  points  along  the  four  curves  defined  by  those 
matrices.    It  then  computes  points  for  each  curve  in  the  U 
direction  using  the  technique  of  forwards  differencing  of  a 
matrix.    Using  these  points  we  can  then  generate  points  along 
the  a  curve  in  the  V  direction. 

void  npatch  (geomx,  geomy,  geomz) 

Coord  geomx[4][4],  geomy  [4]  [4],  geomz  [4]  [4]; 

{ 

register   Point  *  Coord   array [4]; 

/*  One  control  matrix  for  each  curve.  */ 
Matrix  ctrl   ptsl; 
Matrix  ctrl_pts2; 
Matrix  ctrl   pts3; 
Matrix  ctrl   pts4; 
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/*  Load  the  appropriate  control  matrix  for  each  curve.  */ 


/*  Curve  1  */ 


ctrl 

ptsl[0 

[0] 

ctrl 

ptsl[0 

[1] 

ctrl 

ptsl[0 

[2] 

ctrl 

ptsl[l 

[0] 

ctrl 

ptsljl 

[1] 

ctrl 

ptsljl 

[2] 

Ctrl 

ptsl[2 

[0] 

ctrl 

ptsl[2 

[1] 

ctrl 

ptsl[2 

[2] 

ctrl 

ptsl[3 

[0] 

ctrl 

ptsl[3 

[1] 

ctrl 

ptsl[3 

[2] 

/*  Curve  2 

V 

ctrl 

pts2[0 

[0] 

ctrl 

pts2[0 

[i] 

ctrl 

pts2[0 

[2] 

ctrl 

pts2[l 

[0] 

ctrl 

pts2[l 

[1] 

ctrl 

pts2[l 

[2] 

ctrl 

pts2[2 

[0] 

ctrl 

pts2[2 

[1] 

ctrl 

pts2[2 

[2] 

ctrl 

pts2[3 

[0] 

ctrl 

pts2[3 

[1] 

ctrl 

pts2[3 

[2] 

georax 
geomy 
geomz 

georax 
geomy 
geomz 

geomx 
geomy 
geomz 

geomx 
geomy 
geomz 


geomx 
geomy 
geomz 

geomx 
geomy 
geomz 

geomx 
geomy 
geomz 

geomx 
geomy 
geomz 


Oj; 

o]; 

0]; 

0]; 

o]; 

0]: 

o]; 
o]; 

0]: 
01; 


1]: 

i]; 
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/*  Curve  3  */ 

Ctrl   pts3[0][0] 
Ctrl   pts3[0][l] 
Ctrl   pts3[0][2] 

=  geomx[0][2]: 
=  geomy[0][2]; 
=  geomz[0][2]: 

ctrl   pts3[l][0] 
Ctrl   pts3[l][l] 
ctrl   pts3[l][2] 

=  geomx[l][2]; 
=  geomy[l][2]; 
=  geomz[l][2]: 

ctrl   pts3[2][0] 
ctrl   pts3[2][l] 
ctrl   pts3[2][2] 

=  geomx[2][2]; 
=  geomy[2][2]; 
=  geomz[2][2]; 

ctrl   pts3[3][0] 
ctrl   pts3[3][l] 
ctrl   pts3[3][2] 

=  geomx[S][2]: 
=  geomy[3][2]; 
=  geomz[3][2]; 

/*  Curve  4  */ 
ctrl   pts4[0][0] 
ctrl   pts4[0][l] 
ctrl   pts4[0][2] 

=  geomx[0][3]; 

=  geomy[0][3]; 
=  geomz[0][3]: 

ctrl   pts4[l][0] 
Ctrl  pts4[l][l] 
ctrl   pts4[l][2] 

=  geomx[l][3]; 
=  geomy[l][3]; 
=  geomz[l][3]: 

Ctrl  pts4[2][0] 
ctrl   pts4[2][l] 
Ctrl   pts4[2][2] 

=  geomx[2][3]; 
=  geomy[2][3]; 
=  geomz[2][3]; 

ctrl   pts4[3][0] 
ctrl   pts4[3][l] 
Ctrl   pts4[3][2] 

=  geomx[3][3]; 
=  geomy[3][3]; 
=  geomz[3][3]: 
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{ 

register   Point  *  Where: 

register  unsigned  int  Row; 
register  unsigned  int  Column; 
register  unsigned  int  point   count; 
register  unsigned  int  count; 

/*  An  array  of  pointers  to  our  control  matrices.  */ 
float    *matrix  pointer [4]; 

/*  Matrix  used  to  hold  mathematical  results.  */ 
Matrix  interl; 

/*  Initialize  the  array  */ 
matrix   pointer[0]  =  (float  *)  Ctrl   ptsl; 
matrix  pointerfl]  =  (float  *)  ctrl   pts2; 
matrix   pointer[2]  =  (float  *)  ctrl   pts3; 
matrix  pointer[3]  =  (float  *)  ctrl_pts4; 

/*  Get  enough  memory  to  hold  the  points.  */ 
Coord   array [0]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)) 
Coord   arrayjl]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)) 
Coord   array[2]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)) 
Coord   array[3]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)) 

/*  For  each  curve.  */ 
for  (count  =  0;  count  <  4;  countH — \-)  { 

I*  Generate  the  matrix  used  in  the  forward  difference  for  this  curve.  */ 
pushmatrix  (); 

loadmatrix  (matrix  pointer[count]); 
multmatrix  (current   U   basis); 
multmatrix  (Precision   Matrix  U); 
getmatrix  (interl); 
popmatrix  (); 
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/*  For  each  curve  generate  points  in  the  U  parametric  direction.  */ 
for  (point    count  =  0;  point   count  <=  UCURVES;  point    countH — h)  { 

Where  =  (Point  *)  (Coord   array  [count]  +  point   count): 

Where  ->  x  =  interl[3][0 
Where  ->  y  =  interl[3][l 
Where  ->  z  =  interl[3][2] 

/*  Do  the  forward  difference.  */ 
for  (Row  =  3;  Row  >  0;  Row--) 
for  (Column  =  0;  Column  <  4;  Column-(-  +  ) 
inter  1  [Row]  [Column]  =  interl [Row]  [Column]  -f  interl[Row  -  l]  [Column] 

} 


} 

/*  Call  function  to  finish  computations  and  display.  */ 
nspeckle  (Coord    array); 

/*  Return  the  memory  used  back  to  the  system.  */ 
cfree  (Coord   array[0],  (UCURVES  +  1),  sizeof  (Point)) 
cfree  (Coord   array [l],  (UCURVES  +  1),  sizeof  (Point)) 
cfree  (Coord   array [2],  (UCURVES  +  1),  sizeof  (Point)) 
cfree  (Coord   array [3],  (UCURVES  +  1),  sizeof  (Point)) 

} 
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